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ABSTRACT 

The theory for the formation of the first population of stars (Pop III) predicts a IMF composed predominantly 
of high-mass stars, in contrast to the present-day IMF, which tends to yield stars with masses less than 1 M©. 
The leading theory for the transition in the characteristic stellar mass predicts that the cause is the extra cooling 
provided by increasing metallicity and in particular the cooling provided at high densities by dust. The aim 
of this work is to test whether dust cooling can lead to fragmentation and be responsible for this transition. 
To investigate this, we make use of high-resolution hydrodynamic simulations. We follow the thermodynamic 
evolution of the gas by solving the full thermal energy equation, and also track the evolution of the dust 
temperature and the chemical evolution of the gas. We model clouds with diff'erent metallicities, and determine 
the properties of the cloud at the point at which it undergoes gravitational fragmentation. We follow the further 
collapse to scales of an AU when we replace very dense, gravitationally bound, and collapsing regions by a 
simple and nongaseous object, a sink particle. Our results suggest that for metallicities as small as lO'^Z©, dust 
cooling produces low-mass fragments and hence can potentially enable the formation of low mass stars. We 
conclude that dust cooling aff'ects the fragmentation of low-metallicity gas clouds and plays an important role 
in shaping the stellar IMF even at these very low metallicities. 

Subject headings: early universe — hydrodynamics — methods: numerical — stars: formation — stars: lumi- 
nosity function, mass function 



1. INTRODUCTION 

The first burst of star formation in the Universe is thought 
to give rise to massive stars, the so-called Population III, 
with current theory predicting masses in the range 20-150 M© 
dAbel et al .l[2002i Bromm et al., 2002; O'Shea & Norman, 
|2007'i 'Yoshida et al., 2008). This contrasts with present-day 
star formati on, which tends to yield stars with masses less 



than 1 M© ( [Kr ^ 
point in the evolution 



Dupa| [20021 |Chabrier| 
olution of the Univers 



2003 ), and so at some 

Universe there must have been 



a transition from primordial (Pop. Ill) star formation to the 
mode of star formation we see today (Pop. II/I). 
When gas collapses to form stars, gravitational energy is 
transformed to thermal energy and unless this can be dissi- 
pated in some fashion, it will eventually halt the collapse. 
Thermal energy can be dissipated by processes such as atomic 
fine structure line emission, molecular rotational or vibra- 
tional line emission, or the heating of dust grains. In some 
cases, these processes are able to cool the gas significantly 
during the collapse. This te mperature drop can promote grav- 
itational fragmentation (Mac Low & Klessen] |2004[ |Bonnell| 
[et al.^ |200 7 ) by diminishing the Jeans mass, which means 
that instead of forming very massive clumps, with fragment 
masses corresponding to the initial Jeans mass in the cloud, it 
can instead form even more fragments with lower masses. 
If the gas is cooled only by molecular hydrogen emission, nu- 
merical simulations show that the stars should be very mas- 



sive dAbel et al.l [20021 |Bromm et al.l [20021 [O'Shea & No? 



man 



Ho 



2007[ Yoshida et al. 2008). This happens because the 
_ c'ooling becomes inefficient for temperatures bellow 200K 
and densities above lO'^cm"^. At this temperature and density, 
the mean Jeans mass at cloud fragmentation is 1,000 times 
larger than in present-day molecular clouds. 



for an atomic gas with temperature r^ag and number density 

/tfrag • 

Metal line cooling and dust cooling are eff'ective at lower 
temperatures and larger densities, and so the most widely 
accepted cause for the transition from Pop. Ill to Pop. II is 
metal enrichment of the interstellar medium by the previous 
generations of stars. This suggests that there might be a 
critical metallicity Z^nt at which the mode of star formation 
changes. 

The main coolants that have been stu died in the literature 
are CII and 01 fine structure emission fBromm et al.| |2001 



Bromm & Loebl [20031 [Santoro & Shu ll, 2006; Frebe 



et al. 



2007), and dust emission. C and O are identified as the key 



species because in the temperature and density conditions 
that characterise the early phases of Pop. Ill star formation, 
the 01 and CII fine-structure lines dominate over all other 
metal transitions ( jHollenbach & McKee|[1989] ). By equating 
the CII or 01 fine structure cooling rate to the compressional 
heating rate due to free-fall collapse, one can define critical 
abundances [C/H] = -3.5 and [0/H] = -3.(f]for efficient 
metal line cooling (Bromm & Loeb [2003 ). However, previ- 
ous works ( [Jappsen et al.} 2009a, bjshow that this metallicity 
threshold does not represent a critical metallicity: the fact that 
metal-line cooling has a larger value than the compressional 
heating does not necessarily lead to fragmentation. 
Dust-cooling models predict a much lower critical metallicity 
(Zcrit ~ 10"^Z©). The conditions for fragmentation in the 
low-metallicity dust cooling model are predicted to occur in 
high density gas, whe re the di s tances between the fragments 
can be very sm all (Omukai, 2000; Omukai et al., 2005] 
[Schneider etaLl g002, 2006; Schneider & Omukai, 2010). In 
this regime, interactions between fragments will be common. 



rag 



/tfrag 

lO^cm- 



-1/2 



(1) 



1 [X/Y] = logioCNx/Ny)* -logio(Nx/NY)o, for elements X and Y, where 
* denotes the gas in question, and where Nx and Ny are the mass fractions 
of the elements X and Y. 
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and analytic models of fragmentation are unable to predict 
the mass distribution of the fragments. A full 3D treatment, 

following the fragments, is needed. 

Initial attempts we re made by T suribe & Qmukai 



p006| |2008) and jClark et al.| ([200BT However, these 
treatments used a tabulated equation of state, based on results 



from previous one-zone chemical models ( Qmukai et al. 
2QQ5| ), to determine the thermal energy. This approximation 



assumes that the gas temperature adjusts instantaneously to a 
new equilibrium temperature whenever the density changes 
and hence ignores thermal inertia effects. This may yield too 
much fragmentation. 

In this work, we improve upon these previous treatments by 
solving the full thermal energy equation, and calculating the 
dust temperature through the energy equilibrium equation. 
We assume currently that the only significant external heat 
source is the CMB, and include its eff'ects in the calculation 
of the dust temperature. 



2. SIMULATIONS 

2.1. Numerical method 

We model the collapse of a low- metallicity gas c loud using 
a modified version of the Gadget 2 ( [Springel |2005| ) smoothed 
particle hydrodynamics (SPH) code. To enable us to continue 
our simulation beyond the formation of the first very high den- 
sity protoste llar core, we use a sink particle appro ach (Bate, 




1995}, based on the implementation of Jappsen et al. 
k 



Sink particles are created once the SPH particles are 
collapsing, and within an accretion radius, hacc^ which 
is taken to be 1.0 AU. The threshold number density for sink 
particle creation is 5.0 x lO^^cm"^. At the threshold density, 
the Jeans length at the minimum temperature reached by the 
gas is approximately one AU, while at higher densities the 
gas becomes optically thick and begins to heat up. Further 
fragmentation on scales smaller than the sink particle s cale is 
therefore unlikely to occur. For further discussion see Clark 
|etal.| ( |20111 ). 

To treat the chemistry and thermal balance of the gas, we use 
the same approach as in Clark et al. ( 2011| ), with two ad- 
ditions: the inclusion of the effects of dust cooling, as de- 
scribed below, and formation of H2 on the sur face of dust 
grains (see Hollenbach & McKee| [1979]). The IClark et al 



?201 1 ) chemical network and cooling function were designed 
for treating primordial gas and do not include the chemistry 
of metals such as carbon or oxygen, or the effects of cooling 
from these atoms, or molecules containing them such as CO 
or H2O. We justify this approximation by noting t hat previ 



ous studies of very low-metallicity gas (e.g. Q mukai et al. 
[2005 ,2010 ) find that gas-ph ase metals have little influence on 
the thermal state of the gas. Qmukai et al.|p01Q| ) showed that 
H2Q and QH are efficient coolants at 10^ < n < lO^^cm"^ for 
their one-zone model. In their hydrodynamical calculations, 
however, the collapse is faster, and the effect of H2Q and QH 
is not perceptible. Therefore we do not expect oxygen-bearing 
molecules to have a big effect on the thermal evolution of 
the gas. For the metallicities and dust-to-gas ratios consid- 
ered in this study, the dominant sources of cooling are the 
standard primordial coolants (H2 bound-bound emission and 
collision-induced emission) and energy transfer from the gas 
to the dust. 
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Fig. 1 . — Results of our low-resolution simulations, showing the dependence 
of gas and dust temperatures on gas density for metallicities 10~^ and 10~^ 
times the solar value. In red, we show the gas temperature, and in blue the 
dust temperature for the turbulent and rotating cloud. The simple core col- 
lapse is overploted in dark red and green. The points with thinner features 
are from the simulations without rotation or turbulence, while those showing 
more scatter come from the simulations with rotation and turbulence. The 
dashed lines show constant Jeans mass values. 



2.1.1. Dust cooling 

Collisions between gas particles and dust grains can transfer 
energy from the gas to the dust (if the gas temperature T is 
greater than the dust temperature Tgr), or from the dust to the 
gas (if Tgr > T). The rate at which energy is transferred from 
gas to dust is given by ( |Hollenbach & McKee^ |1979j ) 

(2) 



Agr = ^gi-^(TgrVp/(2^r - IkTgr) erg s' 



-1 -3 

cm , 



where rigr is the number density of dust grains, n is the num- 
ber density of hydrogen nuclei, o-gr is the mean dust grain 
cross-section, Vp is the thermal speed of the proton, and / is 
a factor accounting for the ontribution of species other than 
protons, as well as for charge and accommodation effects. We 
assume that o-gr is the same as for Milky Way dust, and that 
the number density of dust grains is a factor Z/Z© smaller than 
the Milky Way value. To compute the rate at which the dust 
grains radiate away energy, we use the approximation (|Sta-| 
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Imatellos et all [2007] ) 



\ A ^ ^cmb^ 
Arad = 4crsb^gr— ^ -p, 



(3) 



where Tcmb is the CMB temperature, (Xsb is the Stefan- 
Boltzmann constant, Kp and kr are the Planck and Rosseland 
mean opacities and S is the column density of gas measured 
along a radial ra y from the particle to the edge of the cloud. 
As explained by Stamatellos et al. (|2QQ7| ), this expression has 
the correct behaviour in the optically thin and optically thick 
limits, and interpolates between these two limits in a smooth 
fashion. In practice, we approximate further by assuming that 
the Planck and Rosseland mean opacities are equaland by us- 
ing the fact that E ~ pLj for a gravitationally collapsing gas, 
where p is the mass density of the gas, and Lj is the Jeans 
length, given by Lj = (ncyGp)^^^, where Cs is the speed of 
sound in the gas. By approximating S in this fashion, we avoid 
the computational difficulties involved with measuring col- 
umn densities directly in the simulation, while still following 
the behaviour of the gas reasonably accurately in the optically 
thick regime. In any case, most of the interesting behaviour 
that we find in our simulations occurs while dust cooling re- 
mains in the optically thin regime. To compute the tempera- 
ture of the dust grains, we assume that the dust is in thermal 
equilibrium, and hence solve the equilibrium equation 

Agr - Arad = 0. (4) 

This equation is transcendental, so we solve it numerically. 

2.1.2. Dust opacity 

We follow the dust opacity model of |Goldsmith| ( |2001| ), 
and we calculate the opacity as a functi on of the du st tem- 
perature in the same fashion as in Ban erjee et al.| |2006). 
To convert from the frequency-dependent opacity given in 
Goldsmith ( 2001| ) to our desired temperature-dependent mean 
opacity, we assume that for dust with temperature Tgr, the 
dominant contribution to the mean opacity comes from fre- 
quencies close to a frequency v that is given by hv = akTg^, 
where a = 2.70. At a reference temperature Tq = 6.75 K, this 
procedure yields an opacity 



^(ro) = 3.3 X 10-26^(n/2pgas) 

= 2.664 X 10-^(1 +4[He]) 



(5) 



where [He] is the helium abundance, and n is the number den- 
sity of hydrogen nuclei. At other temperatures, k oc Tg^., so 
long as Tgr < 200 K. For grain temperatures larger than 200 K, 
it is necessary to account for the eff'ects of ice-mantle evapo- 
ration, while at much higher grain temperatures, the opacity 
falls off' extremely rapidly due to the melting of the grains. 
We account for these eff'ects (see [Semenov et al.| |2003 ) and 
so our opacity varies with dust temperature following the re- 
lationship 

T < 200K 

200K < T < 1500K (6) 
T > 1500K 

2.2. Setup and Initial conditions 



We performed three sets of simulations, two at low resolution 
and one at high resolution. The details are shown in Table [T] 




Resolution 


Number of 


Particle 


Turbulence 


Angular 


Level 


Particles 


Mass 




Momentum 






(10-^Mo) 


(^turb/l^gravl) 


(^rot/l^gravl) 


High 


40 X 10^ 


2.5 


0.1 


0.02 


Low 


4x 10^ 


25.0 


0.1 


0.02 








0.0 


0.00 



TABLE 1 
Simulation properties. 



10^ AU 




Fig. 2. — Number density maps for a slice through the high density region. 
The image shows a sequence of zooms in the density structure in the gas 
immediately before the formation of the first protostar. 

Our low resolution simulations were performed to explore the 
thermal evolution of the gas during the collapse, and had 4 
million SPH particles which was insufficient to fully resolve 
fragmentation. We used these simulations to model the col- 
lapse of an initially uniform gas cloud with an initial num- 
ber density of 10^ cm"^ and an initial temperature of 300 K. 
We modelled two different metallicities (lO'^^Z© and lO'^Z©). 
The initial cloud mass was 1000 M©, and the mass resolu- 
tion was 25 X 10"^ Mq. In one set of low-resolution simula- 
tions the gas was initially at rest, while in the other, we in- 
cluded small amounts of turbulent and rotational energy, with 

^turb/l^gravl = 0.1 and = ^rot/l^gravl = 0.02, whcrC ^grav is 

the gravitational potential energy, ^turb is the turbulent kinetic 
energy and E'rot is the rotational energy. For our high resolu- 
tion simulations, which were designed to investigate whether 
the gas would fragment, we employed 40 million SPH par- 
ticles. We adopted initial conditions similar to those in the 
low-resolution run with turbulence and rotation. Again, we 
simulated two metallicities, lO'^^Z© and lO'^Z©. The mass 
resolution (taken to be 100 times the SPH particle mass) was 
2.5 X IO-^Mq. 

3. ANALYSIS 

3.1. Thermodynamical evolution of gas and dust 

In Figure [TJ we compare the evolution of the dust and gas 
temperatures in the low-resolution simulations. The dust tem- 
perature, shown in the lower part of the panels, varies from the 
CMB temperature in the low density region to the gas temper- 
ature at much higher densities. At densities higher than 10^^- 
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10^^ cm"^, dust cooling starts to be effective and begins to 
cool the gas. The gas temperature decreases to roughly 600 K 
in the 10"^ Z© simulations, and 300 K in the Z = lO'^Z© case. 
This temperature decrease significantly increases the number 
of Jeans masses present in the collapsing region, making the 
gas unstable to fragmentation. The dust and the gas temper- 
atures couple for densities higher then lO^^cm"^, when the 
compressional heating starts to dominate again over the dust 
cooling. The subsequent evolution of the gas is close to adia- 
batic. If we compare the results of the runs with and without 
rotation and turbulence, then the most obvious diff'erence is 
the much greater scatter in the n - T diagram in the former 
case. Variations in the infall velocity lead to diff'erent fluid el- 
ements undergoing difl'erent amounts of compressional heat- 
ing. The overall efl'ect is to reduce both the infall velocity 
and the average compressional heating rate. This allows dust 
cooling to dominate at a density that is up to five times smaller 
than in the case without rotation or turbulence. The gas also 
reaches a lower temperature, cooling down to ^ 200K (instead 
of 300K) for the Z = IQ-^Zq case, and to ^ 400K (instead of 
600K) for the Z = 10"^ Z© case. This behavior shows that it 
is essential to use 3D simulations to follow the evolution of 



the col lapsing gas. A similar eflfect can be seen in ,Clark et al. 
( |2QTT] ). 

3.2. Fragmentation 

We follow the thermodynamical evolution of the gas up to 
very high densities of order lO^^cm"^, where the Jeans mass 
is ^ 10"^M©, and so we need a high resolution simulation to 
study the fragmentation behaviour. The transport of angular 
momentum to smaller scales during the collapse leads to the 
formation of a dense disk-like structure, supported by rotation 
which then fragments into several objects. Figure |2] shows the 
density structure in the gas immediately before the formation 
of the first protostar. The top-left panel shows a density slice 
on a scale comparable to the size of the initial gas distribu- 
tion. The structure is very filamentary and there are two main 
overdense clumps in the center. If we zoom in on one of the 
clumps, we see that its internal structure is also filamentary. 
We can follow the collapse down to scales of the order of an 
AU, but at this point we reach the limit of our computational 
approach: as the gas collapses further, the Courant timestep 
becomes very small, making it difficult to follow the further 
evolution of the cloud. In order to avoid this difficulty, we 
replace very dense, gravitationally bound, and collapsing re- 
gions by sink particles. Once the conditions for sink particle 
creation are met, they start to form in the highest density re- 
gions (Figure [3]). Due to interactions with other sink particles 
that result in an increase in velocity, some sink particles can 
be ejected from the high-density region, but most of the par- 
ticles still remain within the dense gas. Within 137 years of 
the formation of the first sink particle, 45 sink particles have 
formed. At this time, approximately 4.6M© of gas has been 
accreted by the sink particles. 

3.3. Properties of the fragments 

Figure H] shows the mass distribution of sink particles when 
we stop the calculation. We typically find masses below IM©, 
with somewhat smaller values in the 10""^ Z© case compared 
to the 10"^Z© case. Both histograms have the lowest sink par- 
ticle mass well above the resolution limit of 0.0025M©. Note 
that in both cases, we are still looking at the very early stages 
of star cluster evolution. As a consequence, the sink particle 
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Fig. 3. — Number density map showing a slice in the densest clump, and the 
sink formation time evolution, for the 40 million particles simulation, and Z 
= 10~^Zo. The box is lOOAU x lOOAU and the time is measured from the 
formation of the first sink particle. 

masses in Figure |4] are not the same as the final protostellar 
masses - there are many mechanisms that will aff'ect the mass 
function, such as continuing accretion, mergers between the 
newly formed protostars, feedback from winds, jets and lu- 
minosity accretion, etc. Nevertheless, we can speculate that 
the typical stellar mass is similar to what is observed for Pop 
II stars in the Milky Way. This suggests that the transition 
from high-mass primordial stars to Population II stars with 
mass function similar to that at the present day occurs early 
in the metal evolution history of the universe, at metallicities 
Zcrit < 10"^Z©. The number of protostars formed by the end 
of the simulation, for both 10~^ Zq (45) and 10"^ Z© (19) cases, 
is much larger than the initial number of Jeans masses (3) in 
the cloud. 
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Fig. 4. — Sink particle mass function at the end of the simulations. High 
and low resolution results and corresponding resolution limits are shown. To 
resolve the fragmentation, the mass resolution should be smaller than the 
Jeans mass at the point in the temperature-density diagram where dust and 
gas couple and the compressional heating starts to dominate over the dust 
cooling. At the time shown, around 5 M© of gas had been accreted by the 
sink particles in each simulation. 

4. CONCLUSIONS 

In this paper we have addressed the question of whether 
dust cooling can lead to the fragmentation of low-metallicity 
star-forming clouds. For this purpose we performed numeri- 
cal simulations to follow the thermodynamical and chemical 
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evolution of collapsing clouds. The chemical model included 
a primordial chemical network together with a description 
of dust evolution, where the dust temperature was calculated 
by solving self-consistently the thermal energy equilibrium 
equation. 

We performed three sets of simulations, two at low resolution 
and one at high resolution (Table [T]). All simulations had 
an initial cloud mass of 1000 M©, number density of 10^ 
cm"^, and temperature of 300K. We tested two different 
metallicities (lO'^^Z© and lO'^Z©), and also the inclusion of 
small amounts of turbulent and rotational energies. 
We found in all simulations that dust can effectively cool 
the gas, for number densities higher than 10^ ^cm"^. An 
increase in metallicity implies a higher dust-to-gas ratio, and 
consequently stronger cooling by dust. This is reflected in a 
lower temperature of the dense gas in the higher metallicity 
simulation. 

For the low resolution case, we tested the effect of adding 
turbulence and rotation. These diminish the infall velocity, 
leading to different fluid elements undergoing different 
amounts of compressional heating. This lack of heating 
allows the gas to reach a lower temperature. 
We found that the transport of angular momentum to smaller 
scales lead to the formation of a disk-like structure, which 
then fragmented into a number of low mass objects. 



We conclude that the dust is already an efficient coolant 
even at metallicities as low as 10"^ or IO'^^Zq, in agreement 
with previous works ( Clark et al. , "2008 ; "Qmukai et al. , 2010} 
Schneider et al. , 2002 , 2006 ; Tsuribe & Omukai , 2006 , 20^ 
Our results support the idea that dust cooling can play an 
important role in the fragmentation of molecular clouds and 
the evolution of the stellar IMF. 
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